\chapter{Model Equations and Discretization} \label{Chapter2}

\section{The Variational Boussinesq Model (VBM)}
The basic idea of the model is to minimize the pressure in the whole fluid. Starting point of the derivation of the model equation comes from the physical model of the Navier-Stokes Equations. The Euler equations are a special case of the Navier Stokes equations where viscous forces are neglected. For ocean waves, this is considered to be a valid assumption as the Reynolds number is typically large, and thus surface waves in water are a superb example of a stationary and ergodic random process. Also, rotational effects are negligible in large water waves.
The Euler equation is given by: 

\begin{equation} \label{EulerEquation}
\frac{\partial \textbf{u}}{\partial t} + (\textbf{u}.\nabla)\textbf{u} = \frac{1}{\rho}\nabla p -\textbf{g},
\end{equation}

The fluid velocity is given by \textbf{u}, $\rho$ denotes the mass density, $p$ the pressure and \textbf{g} the gravitational constant.
Understanding  physics behind the problem is important as it provides information on the numerical techniques required for the solution of the problem.
Note that the inviscid and irrotational assumption is not necessarily valid near solid boundaries, where very small flow structures associated with turbulence result from the no-slip boundary condition. In the current model, this has been ignored and the assumptions are considered to be valid throughout the domain.

Equation $\ref{EulerEquation}$ is converted to the instationary Bernoulli equation (fluid being assumed irrotational) , and then the pressure is integrated over the whole domain. The basic idea of the variational Boussinesq model is to minimize the total pressure  $p$. For irrotational flows, velocity can be represented by a velcoity potential ($\phi$). Another parameter $\zeta$ (water level) is also introduced while performing integration. The problem thus becomes finding these functions ($\phi$ and $\zeta$). The vertical structure of the flow is often known. The velocity potential can be written as a series expansion in predefined vertical shape functions, thus reducing the 3D model to a 2D model.

\subsection{Hamiltonian Description}
Being an irrotational and incompressible flow, the system of equations has a Hamiltionian structure, where the force acting on the system is the gravitational and the pressure force. The Hamiltionian $H(\zeta,\varphi)$ is given by the sum of kinetic and potential Energy. In order to obtain the minimal value of total pressure, the Hamiltionian function is minimized with respect to variables $\zeta$ and $\varphi$ \cite{Wout}. This results in :

\begin{subequations}
\label{Hamiltonian}
\begin{equation}
 \dfrac{\partial \zeta}{\partial t} = \nabla H_{\varphi}(\varphi, \zeta)
 \end{equation}
 \begin{equation}
 \dfrac{\partial \varphi}{\partial t} = \nabla H_{\zeta}(\varphi, \zeta)
\end{equation}
\end{subequations}

where $\nabla H_{\zeta}$ refers to the partial derivative of $H$ with respect to $\zeta$.
The structure of the above equations will play an important role in the discussion of temporal integration which is discussed in Chapter 4.

\subsection{Linearized Model Description}
The variational Boussinesq model constructed based on the above simplifications is non-linear in nature. In order to simplify these equations and thus reducing the computational effort, the Hamiltionian equations given by \Cref{Hamiltonian} is linearized. Detailed derivation from the Hamiltionian to linear set of equations is described in \cite{Wout}.

The final set of equations after performing the linearization are given by:

\begin{subequations}
\label{basicEqns}
\begin{equation}
\label{2.2a}
\frac{\partial \zeta}{\partial t} + \nabla .(\zeta\textbf{U}+h\nabla \varphi - h\mathcal{D}_{0}\nabla \psi) = 0,
\end{equation}
\begin{equation}
\label{2.2b}
\frac{\partial \varphi}{\partial t} + \textbf{U}.\nabla \varphi + g\zeta = P_{s},
\end{equation}
\begin{equation}
\label{2.2c}
\mathcal{M}_{0}\psi + \nabla .(h\mathcal{D}_{0}\nabla \varphi - \mathcal{N}_{0}\nabla \psi) = 0,
\end{equation}
\end{subequations}

Notes:
\begin{enumerate}
 \item Equations \ref{basicEqns} are solved for three basic variables: water level $\zeta$, surface velocity potential (2-D) $\varphi$ and vertical structure $\psi$. Splitting of the velocity potential into the surface potential and the vertical shape function is given by:

\begin{equation}
\phi (x,y,z,t) = \varphi (x,y,t) + f(z)\psi (x,y,t)
\end{equation}

\item Shape function $f$ is chosen among either a parabolic or a cosine-hyperbolic shape. The model parameters (functionals $\mathcal{D, M}$ and $\mathcal{N}$) are computed using the shape function $f$.
\item The water depth $h =h(x,y,t)$ is relative to the reference level $z=0$. The bottom of a basin, river or sea is thus at level $z=-h$.
\item The total velocity can be split into the average current velocity and the velocity due to the wave front. $\textbf{U} = \textbf{U}(x,y,t)$ is the time average horizontal velocity of the current and is used as an external input in the model.
\item Equations \ref{basicEqns} represent the motion of waves. The impact of a moving ship is not seen directly. In order to model a ship, a pressure pulse on the water surface is defined. In the pressure functional given by \Cref{2.2b}, $P_{s}$  represented the source term with $P_{s} := -\frac{p_{s}}{\rho}$. The draft of a ship's hull is the vertical distance between the waterline and the bottom of the hull. Given the draft of a ship, $p_{s}$ is computed as the hydrostatic pressure at the given depth. Let the draft be given as $d_s$, then $p_s = g d_s \alpha (x,y)$ with $\alpha(x,y)$ a shape function with one in the middle of the ship and zero on the horizontal boundary of the ship. Alternatively, the shape of the ship's hull can be imported from an external surface description file.
\end{enumerate}

\section{Numerical Discretization}

\subsection{Computational Domain}
Both Elwin and Martijn have assumed the domain to be rectangular, and divided the domain in a rectilinear Cartesian grid. The dimensions of the domain are $L_{x}\times L_{y}$. It is divided into $N_{x}\times N_{y}$ grid points. The outermost nodes represent the physical boundary. The mesh spacing in each direction is given as $\Delta x = \frac{L_{x}}{N_{x}-1}$ and $\Delta y = \frac{L_{y}}{N_{y}-1}$. An example is given in  \ref{figure2.1}

\begin{figure}[ht]
\label{figure2.1}
\centering
\includegraphics[width=1.0\textwidth]{grid}~\\[1cm]
\caption{Physical domain and corresponding Cartesian grid}
\end{figure}

\subsection{Spatial Discretization}
\label{section_spatial}
The model Equations \ref{basicEqns} are discretized using a finite volume method. The variables are evaluated at the grid points and the finite volumes are rectangles of size $\Delta x \times \Delta y$ centered around the grid point. The derivatives are approximated with centered differences yielding a five-point stencil. For the grid point located at C (= center) the surrounding control volume V and its four nearest neighbors (N  = north,E = east, S = south, W = west) are indicated in \Cref{figure2.1}. On node $(i,j)$, the discretized versions of the variables $\zeta , \varphi \text{ and } \psi $ are given by $\zeta _{ij} , \varphi _{ij}  \text{ and } \psi _{ij}$. In order to put the variables in matrix format, one dimensional ordering of the variables is defined which gives the vector $\vec{\zeta}, \vec{\varphi} \text{ and }  \vec{\psi}$. The spatial discretization can be written as :

\begin{equation}
\frac{\mathrm{d}}{\mathrm{dt}} \begin{bmatrix} \vec{\zeta} \\ \vec{\varphi} \\0 \end{bmatrix} +
\begin{bmatrix} S_{\zeta \zeta} & S_{\zeta \varphi} & S_{\zeta \psi} \\ S_{\varphi \zeta} & S_{\varphi \varphi} & S_{\varphi \psi} \\S_{\psi \zeta} & S_{\psi \varphi} & S_{\psi \psi}  \end{bmatrix} \begin{bmatrix} \vec{\zeta} \\ \vec{\varphi} \\ \vec{\psi}  \end{bmatrix} = \begin{bmatrix} 0 \\ P_s \\ 0 \end{bmatrix} 
\end{equation}

The matrix S's are given by five point stencils as follows:

\begin{equation}
\label{Stencil_zz}
 S_{\zeta \zeta} : \begin{bmatrix} 0 && \frac{1}{2 \Delta y}\overline{V_N} && 0\\
 -\frac{1}{2 \Delta x}\overline{U_W}    &&  \frac{1}{2 \Delta y}\overline{V_N} -\frac{1}{2 \Delta y}\overline{V_S}   -\frac{1}{2 \Delta x}\overline{U_W} +\frac{1}{2 \Delta x}\overline{U_E} &&   \frac{1}{2 \Delta x}\overline{U_E} \\
 0 && -\frac{1}{2 \Delta y}\overline{V_S} && 0  \end{bmatrix} 
\end{equation}

where $U$ and $V$ denotes the current velocity in $x$ and $y$ direction respectively.


\begin{equation}
 S_{\zeta \varphi} : \begin{bmatrix} 0 && -\frac{1}{\Delta y^2}\overline{h_N} && 0\\
 -\frac{1}{\Delta x^2}\overline{h_W}    && \frac{1}{\Delta y^2}\overline{h_N}+ \frac{1}{\Delta x^2}\overline{h_W}+ \frac{1}{\Delta x^2}\overline{h_E}+ \frac{1}{\Delta y^2}\overline{h_S} &&   -\frac{1}{\Delta x^2}\overline{h_E} \\
 0 && -\frac{1}{\Delta y^2}\overline{h_S} && 0  \end{bmatrix} 
\end{equation}


\begin{equation}
 S_{\zeta \psi} : \begin{bmatrix} 0 && -\frac{1}{\Delta y^2}\overline{h_N}\overline{D_N} && 0\\
 -\frac{1}{\Delta x^2}\overline{h_W}\overline{D_W}    && \frac{1}{\Delta y^2}( \overline{h_N}\overline{D_N}+\overline{h_S}\overline{D_S} )+ \frac{1}{\Delta x^2}( \overline{h_W}\overline{D_W} + \frac{1}{\Delta x^2}\overline{h_E}\overline{D_E}) && -\frac{1}{\Delta x^2}\overline{h_E} \overline{D_E} \\
 0 && -\frac{1}{\Delta y^2}\overline{h_S} \overline{D_S} && 0  \end{bmatrix} 
\end{equation}


\begin{equation}
 S_{\varphi \zeta} = g , \quad S_{\varphi \psi}  = 0, \quad S_{\psi \zeta} = 0
\end{equation}


\begin{flalign}
 & S_{\varphi \varphi} : \begin{bmatrix} 0 && \frac{1}{2 \Delta y}\overline{V_N} && 0\\
 -\frac{1}{2 \Delta x}\overline{U_W}    &&  -(\frac{1}{2 \Delta y}\overline{V_N} -\frac{1}{2 \Delta y}\overline{V_S} -\frac{1}{2 \Delta x}\overline{U_W} +\frac{1}{2 \Delta x}\overline{U_E}) &&   \frac{1}{2 \Delta x}\overline{U_E} \\
 0 && -\frac{1}{2 \Delta y}\overline{V_S} && 0  \end{bmatrix} 
\end{flalign}
 
 
 \begin{equation}
 S_{\psi\varphi} : \Delta x \Delta y \begin{bmatrix} 0 && \frac{1}{\Delta y^2}\overline{h_N}\overline{D_N} && 0\\
 \frac{1}{\Delta x^2}\overline{h_W}\overline{D_W}    && -(\frac{1}{\Delta y^2}\overline{h_N}\overline{D_N}+ \frac{1}{\Delta x^2}\overline{h_W}\overline{D_W} + \frac{1}{\Delta x^2}\overline{h_E}\overline{D_E} +\frac{1}{\Delta y^2}\overline{h_S}\overline{D_S}) && \frac{1}{\Delta x^2}\overline{h_E} \overline{D_E} \\
 0 && \frac{1}{\Delta y^2}\overline{h_S} \overline{D_S} && 0  \end{bmatrix} 
\end{equation}


\begin{equation}
\label{S_varphi}
 S_{\psi \psi} : \Delta x \Delta y \begin{bmatrix} 0 && -\frac{1}{\Delta y^2}\overline{\NN_N} && 0\\
 -\frac{1}{\Delta x^2}\overline{\NN_W}    && \frac{1}{\Delta y^2}\overline{\NN_N}+ \frac{1}{\Delta x^2}\overline{\NN_W}+ \frac{1}{\Delta x^2}\overline{\NN_E}+ \frac{1}{\Delta y^2}\overline{\NN_S} + \mathcal{M} &&   -\frac{1}{\Delta x^2}\overline{\NN_E} \\
 0 && -\frac{1}{\Delta y^2}\overline{\NN_S} && 0  \end{bmatrix} 
\end{equation}
 

The system can be written as :
\begin{subequations}
 \begin{equation}\label{timeeqn} 
  \dot{\textbf{q}} = L\textbf{q}+f,
 \end{equation}
 \begin{equation}\label{spatialeqn}  
  S\vec{\psi} = \textbf{b}
 \end{equation}
\end{subequations}

with $\bq= \begin{bmatrix} \vec{\zeta} \\ \vec{\varphi} \end{bmatrix}$ and $\dot{\bq}$ its time derivative. 
The matrix $L = -\begin{bmatrix} S_{\zeta \zeta} & S_{\zeta \varphi} \\ S_{\varphi \zeta} & S_{\varphi \varphi}\end{bmatrix}$ is the spatial discretization matrix and $\textbf{f} = - \begin{bmatrix}  S_{\zeta \psi}\vec{\psi} \\ S_{\varphi \psi}\vec{\psi} \end{bmatrix} $.

Elwin and Martijn focused on solving the system of equations represented by \Cref{spatialeqn}. One of the tasks of the current literature review is to study various time integration techniques for \Cref{timeeqn}.

\subsection{Temporal Integration :Leap Frog Scheme }
\label{section_leap_frog}
In the current work by Elwin and Martijn, the leapfrog method has been used to integrate equation \ref{timeeqn}. The Leapfrog method is one of the Symplectic Integration techniques designed for the numerical solution of Hamilton's equations given by \Cref{Hamiltonian}. More about Symplectic Integration is discussed in Chapter 4.

The method described in \cite{Wout} is an explicit scheme and depends on two previous time steps (a so-called multistep) method. The exact solution to \Cref{timeeqn} is approximated at the time intervals $t_n=n\Delta t$, $n=0,1,2,\dots $ with $\Delta t >0$ being the time step size. The numerical approximations are denoted by $\bq^n \approx \bq (t_n)$.

To keep the derivation short, we will first focus on the fixed constant step size $\Delta t := t_{n+1} - t_n$. A Taylor series expansion gives:

\begin{subequations}\label{TaylorSeries} 
 \begin{equation}
  \bq^{n+1} = \bq^n +  \Delta t\dot{\bq}^n + \dfrac{1}{2}\Delta t^2 \ddot{\bq}^n + \dfrac{1}{6}\Delta t^3 \dddot{\bq}^n + \mathcal{O}(\Delta t^4) ,
 \end{equation}
 \begin{equation}
  \bq^{n-1} = \bq^n -  \Delta t\dot{\bq}^n + \dfrac{1}{2}\Delta t^2 \ddot{\bq}^n - \dfrac{1}{6}\Delta t^3 \dddot{\bq}^n + \mathcal{O}(\Delta t^4)
 \end{equation}
\end{subequations}

Subtracting the second equation from first, we obtain:

 \begin{equation}
  \bq^{n+1} -\bq^{n-1}  = 2 \Delta t\dot{\bq}^n + \mathcal{O}(\Delta t^3) ,
 \end{equation}
 
which reduces to 
 
  \begin{equation}
  \dot{\bq}^n  = \dfrac{\bq^{n+1} -\bq^{n-1}}{2 \Delta t} + \mathcal{O}(\Delta t^2) ,
 \end{equation}
 
The explicit nature limits the size of the time-step for stability reasons. For example, given the equation $\dot{y} = \lambda y $ ($ \lambda $ arises as an eigenvalue of a local Jacobian and could be complex), the leapfrog method is stable only for $|\lambda \Delta t \leq 1|$. In order to approximate the interaction between the ship and the waves, a grid size of the order of half a meter or smaller is required. This limits the time step to $0.01$ seconds maximum, while the maneuvering simulator at MARIN often can use time steps as large as $0.1$ seconds. This provides the motivation to explore other time integration methods like the implicit methods, which would allow us to use larger time steps, without causing instability.
